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Abstract 


I 

\f 

'  '■  The  growth  rate  of  the  hydromagnetic  Rayleigh-Taylor 
instability  is  approximated  here  for  an  accelerating  plasma 
slab.  The  slab  is  chosen  as  a  large-radius  approximation 
to  an  imploding  cylindrical  foil.  A  normal  mode  solution 
of  the  MHD  equations  is  assumed,  resulting  in  an  integral 
relation  for  the  instability  growth  rate.  The  Rayleigh- 
Ritz  variational  method  is  applied  to  the  relation  to 
.  estimate  the  growth  rate.  A  linearly  decreasing  magnetic 
field  is  assumed  in  the  slab  perpendicular  to  the  accel¬ 
eration.  A  corresponding  equilibrium  mass  density  profile 
is  then  found.  Growth  rate  estimates  are  then  made  for 
these  profiles.  Calculations  are  made  for  perturbation 
wavevectors  perpendicular  to  the  acceleration  and  at  an 
angle  e  to  the  magnetic  field.  The  growth  rates  for 
B  =  90^  compare  favorably  with  LeLevier  et  al's  results 
for  a  continuous  density  transition  (14)  .  Growth  rates 
for  8  =  0*  are  stable  for  all  perturbation  wavelengths 
and  magnetic  field  strengths.  This  contradicts  prior 
results  in  both  slab  and  cylindrical  geometry  and  suggests 
an  error  in  this  work, 
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I.  Introduction 


In  this  thesis  the  growth  rates  of  the  hydromagnetic 
Rayleigh-Taylor  instability  for  an  accelerating  plasma  slab 
are  approximated  using  the  Rayleigh-Ritz  method.  An  accel¬ 
erating  slab  is  chosen  as  an  approximation  to  an  imploding 
cylindrical  shell.  Discussion  of  a  current  experimental 
effort  points  out  the  usefulness  of  this  work. 


The  Air  Force  Weapons  Laboratory  is  currently  generat¬ 
ing  short,  intense  x-ray  pulses  with  the  SHIVA  machine. 
SHIVA  consists  of  a  megajoule  capacitor  bank  connected  to 
a  pair  of  electrodes.  A  thin  (about  one  micron  thick) 
cylindrical  foil  is  placed  between  the  electrodes.  This 


•» 

-*  * 


configuration  is  schematically  illustrated  in  Figure  1. 

When  the  capacitor  bank  is  switched  across  the  electrodes 
the  foil  ionizes  and  expands  into  a  plasma  sheath.  The 
mega-amp  currents  passing  through  the  foil  generate  a  mega¬ 
gauss  azimuthal  magnetic  field  outside  of  the  foil.  This 
external  field  is  equivalent  to  that  generated  by  a  wire 
placed  on  axis : 


B 


(1.1) 


There  is  now  a  5  x  3  force  directed  radially  inward  which 
causes  the  sheath  to  implode  towards  the  axis .  With  the 
large  current  and  magnetic  field,  the  sheath  obtains  an 
average  velocity  on  the  order  of  107  cm/sec  as  it  collapses 
on  axis.  The  sheath  kinetic  energy  is  then  thermalized, 
producing  a  hot,  dense  plasma.  This  plasma  radiates 
strongly  in  the  x-ray  region  (see  Figure  1) .  X-ray  powers 
on  the  order  of  10 ^ 2  watts  have  already  been  observed. 

Experience  shows,  however,  that  SHIVA  performance  in 
certain  configurations  is  degraded  by  the  presence  of 
magnetohydrodynamic  (MHD)  instabilities.  To  see  this, 
transform  to  the  rest  frame  of  the  plasma  sheath's  center 
of  mass.  The  plasma  is  then  being  accelerated  against  a 
magnetic  pressure 


P 


m 


(1.2) 


This  is  an  unstable  situation  known  as  the  hydromagnetic 
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Rayleigh-Taylor  instability  (13) .  This  instability  causes 
characteristic  Rayleigh-Taylor  ripples  to  grow  along  the 
outer  plasma  boundary.  These  ripples  grow,  and  numerical 
simulations  show  that  they  can  destroy  the  sheath  if  the 
wavelength  of  the  ripples  approaches  the  sheath  thickness 
(17:274).  This  breakup  in  turn  hinders  the  plasma  thermal- 
ization  by  increasing  the  sheath  thickness  and  therefore 
reduces  the  x-ray  output.  A  good  understanding  of  the 
Rayleigh-Taylor  growth  rates  for  the  SHIVA  implosion  is 
therefore  important. 

Scope 

This  work  is  intended  to  approximate  the  Rayleigh- 
Taylor  growth  rates  for  the  early  part  of  the  implosion, 
when  the  linear  Rayleigh-Taylor  instability  dominates. 

Later  in  the  implosion,  the  instability  evolves  into  a  non¬ 
linear  stage  which  is  more  amenable  to  computer  simulation 
(9) .  Growth  rates  for  azimuthal  perturbations  are  espe¬ 
cially  interesting  since  the  current  computer  simulations 
operate  in  an  r-z  geometry  only.  The  perturbations  con¬ 
sidered  will  have  both  z  and  0  variation.  The  effect  of 
the  azimuthal  magnetic  field  as  it  diffuses  into  the  plasma 
sheath  will  also  be  considered.  This  field  acts  to  help 
stabilize  the  instability.  This  problem  will  be  attacked 
here  in  a  planar  approximation  with  approximate  forms  for 
the  magnetic  field  and  plasma  density  profiles.  These 
9*  approximations  should  be  good  in  the  region  of  interest 
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near  the  outer  edge  of  the  sheath  (where  the  instability 
occurs) .  Finally,  these  growth  rates  will  be  approximated 
using  a  variational  technique  (the  Rayleiah-Ritz  method) 
because  an  exact  solution  could  not  be  found.  An  exact 
solution  for  the  limiting  case  of  z*-axis  perturbations  has 
recently  been  reported  by  Tsai,  Liskow,  and  Wilcox  (19) 
for  essentially  the  same  planar  problem.  They  also  were 
forced  to  use  numerical  methods  to  handle  azimuthal  pertur¬ 
bations  in  the  sheath,  however. 

Assumptions 

All  work  here  was  done  in  planar  geometry  with  a  thin 
plasma  sheath.  Although  the  actual  geometry  is  cylindrical, 
the  agreement  between  the  two  geometries  should  be  good 
during  the  early  stages  of  the  implosion.  Cylindrical 
convergence  effects  are  not  yet  significant  (10) ,  making 
a  planar  approximation  acceptable. 

Several  assumptions  are  made  about  the  sheath  plasma 
itself.  The  sheath  is  assumed  to  be  incompressible,  to  be 
describable  by  the  ideal  MHD  equations,  and  to  have  zero 
electrical  resistivity,  all  essentially  for  convenience. 

It  should  be  mentioned,  however,  that  a  magnetic  field 
cannot  diffuse  into  an  infinitely  conducting  plasma.  Mag¬ 
netic  fields  will  therefore  be  frozen  into  the  plasma  sheath 


under  the  assumption  that  the  fields  would  physically  have 
diffused  into  the  plasma  and  then  have  been  frozen  into 
place. 


The  next  assumption  is  that  the  plasma  viscosity  is 
zero.  Turchi  and  Baker  state  that  for  a  typical  SHIVA  con¬ 
figuration  the  viscous  Reynolds  number  exceeds  10$.  The 
size  of  this  quantity  implies  that  viscosity  can  be  ignored 
inside  the  sheath  (20:4943). 

Organization 

The  body  of  this  work  is  organized  as  follows :  In 
Chapter  II  the  MHD  equations  are  presented.  These  equations 
are  linearized  and  a  normal  mode  solution  is  assumed, 
resulting  in  an  integral  relation  for  the  growth  rates  of 
the  Rayleigh-Taylor  instability.  This  relation  cannot  be 
solved  exactly,  so  the  Rayleigh-Ritz  variational  principle 
is  presented  as  an  approximate  method  of  solution.  Approx¬ 
imate  forms  for  the  magnetic  field  and  plasma  density  pro¬ 
files  in  the  plasma  sheath  are  developed  in  Chapter  III  for 
insertion  into  the  integral  relation.  The  Rayleigh-Ritz 
method  is  then  applied  to  the  integral  relation  to  approx¬ 
imate  the  growth  rates.  To  do  this,  a  computer  program 
was  written.  Its  results  are  reported  in  Chapter  IV. 

This  work  is  done  in  CGS-EMU  units. 
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II.  Background  Theory 

Application  of  the  MHD  Equations  (2:48-56,  3:428-66) 

The  magnetohydrodynamic  (MHD)  model  is  one  of  the 

simplest  models  of  a  plasma  available.  It  considers  the 

plasma  to  be  a  single  conducting  fluid.  In  an  ideal  MHD 

model,  the  fluid  is  assumed  to  be  infinitely  conducting. 

This  in  turn  implies  local  charge  neutrality  in  the  plasma. 

It  is  also  assumed  here  that  the  fluid  is  incompressible 

and  has  no  material  viscosity. 

% 

We  wish  here  to  model  a  plasma  foil  with  a  center  of 
mass  acceleration  g.  It  is  therefore  convenient  to  work 
in  an  accelerating  reference  frame  moving  with  the  plasma 
center  of  mass.  A  uniformly  accelerating  frame,  the  MHD 
equations,  are  identical  to  those  in  an  inertial  coordinate 
system  except  for  the  addition  of  a  +  pg  force  term  to  the 
equation  of  motion  (19:1677).  If  the  center  of  mass  veloc¬ 
ity  is  much  less  than  the  speed  of  light,  which  is  certainly 
the  case  for  an  accelerating  plasma,  the  current  density 
J  and  magnetic  field  B  have  essentially  the  same  magnitude 
in  both  the  laboratory  and  plasma  frames.  The  necessary  MHD 
equations  are  then  as  follows: 

1.  The  equation  of  motion 

p  =  ~  7Pmat  +  c  ^  x  §  +  pg  (2.1) 

2.  The  incompressibility  equation 

9  •  V  =  0  (2.2) 
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Figure  2.  (a)  The  cylindrical  implosion 
geometry  in  the  lab  frame,  (b)  The 
planar- geometry  under  consideration  in 
the  rest  frame  of  the  plasma  center  of 

mass. 


3.  The  mass  continuity  equation 


|f  +  $  •  Vp  =  0  (2.3) 

and  4.  Ampere's  Law  (with  the  displacement  current 
ignored) 

V  x  3  =  5  (2.4) 

c 

By  inserting  (2.4)  into  (2.1),  using  a  vector  identity  for 
(v  x  S)  x  S  and  assuming  that  S  J_  V§,  the  Lorentz  force  can 
be  replaced  by 

f  3  x  §  =  -  +  -j-  §  •  v§ 

c  8ir  4ir 

=  -VPmag  +  £7  S  •  vS  (2.5) 

The  Lorentz  force  therefore  consists  here  of  a  pres¬ 
sure  gradient  component  perpendicular  to  the  magnetic  field 
lines  and  a  restoring  force  component  opposing  bending  of 
the  magnetic  field  lines.  If  the  pressure  Pmat  i-n  (2.1) 
is  replaced  by  a  total  pressure  P  =  Pmat  +  pmag»  (2.1) 
becomes 

=  -  VP  +  pg  +  S  •  vS  (2.6) 

These  equations  are  here  applied  to  the  geometry  shown 
in  Figure  2  with  p  =  p(z),  P  =  P(z),  and  g  =  -gz.  The 
acceleration  is,  of  course,  not  constant  in  a  cylindrical 
implosion.  The  assumption  here  of  a  constant  acceleration 
is  equivalent  to  examining  the  sheath  at  one  particular 
instant  during  the  implosion.  Acceleration  in  the  +z 
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direction  in  Figure  2b  corresponds  to  an  acceleration  in 
cylindrical  coordinates  in  the  -r  direction  (see  Figure 
2a) .  Figure  2b  therefore  corresponds  to  an  implosino  in 
the  positive  z  direction. 

The  equations  (2.2),  (2.3),  and  (2.6)  are  first 
linearized,  with  variables  <J>  written  as  the  sum  of  an 
equilibrium  value  <f>0  and  a  small  perturbed  value  <f>^  which 
contains  all  the  time  dependence.  The  density  pressure 
and  velocity  are  therefore  written  as  p  =  pQ  +  p^,  P  =  PQ 
£  P^,  and  ^  +  V^.  respectively.  From  this  point  on 

the  equilibrium  values  will  be  written  without  subscripts. 
When  these  forms  for  the  variables  are  inserted  into  the 
equations,  quadratic  terms  in  perturbed  quantities  are 
ignored  and  time  derivatives  of  equilibrium  values  equal 
zero.  The  equilibrium  velocity  V0  =  0  in  the  frame  of  the 
sheath. 

A  normal  mode  assumption  is  now  made.  All  perturbed 
quantities  are  assumed  to  have  the  form 

$1  =  (some  function  of  z) .  exp  (ikxx  +  ikvy  +  ft) 

Y  (2.7) 

dependence,  with  k2  =  kx2  +  ky2.  Any  arbitrary  perturba¬ 
tion  can  then  be  represented  as  a  Fourier  series  of  expo¬ 
nential  terms.  Only  the  fastest  growing  Fourier  component 
is  important,  however.  The  simple  dependence  of  (2.7)  is 
therefore  adequate  here.  Inserting  the  dependence  in  (2.7) 
into  the  linearized  equations  implies  the  substitutions 
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(2.8) 


d  . ,  d  . ,  ,  d 

dx  lkx'  dy  ■"  lkY'  and  dt  "*■  Y 

The  only  remaining  derivatives  are  with  respect  to  z. 

The  details  of  the  linearization  and  normal  mode  solu¬ 
tion  are  carried  out  by  Chandrasekhar  (3:428-430,  457-466) 
for  both  the  hydrodynamic  and  hydromagnetic  cases  (the  non¬ 
conducting  and  conducting  cases,  respectively).  In  the 
hydrodynamic  case  the  result  is  the  differential  equation 

D(pDw)  -  pk2w  +  ^2  (Dp)  w  =  o  (2.9) 

» 

where  D  =  and  w  =  V2 .  This  equation  is  in  Sturm- 
Liouville  form.  Sturm-Liouville  theory  states  (3:432)  that 
if  Dp  is  positive  everywhere,  then  all  eigenvalues  y2  are 
positive.  If  Dp  is  everywhere  negative,  all  eigenvalues 
are  negative.  Finally,  if  Dp  is  anywhere  positive  then 
there  exists  at  least  one  positive  eigenvalue.  The  time 
dependence  (2.7)  means  that  an  initial  perturbation  will 
grow  exponentially  if  y2  >  0,  resulting  in  an  unstable 
state.  If  y  <  0,  the  amplitude  of  a  perturbation  oscil¬ 
lates  periodically,  yielding  a  stable  state.  Combining 
these  facts,  we  see  that  the  configuration  is  unstable  if 
Dp  is  positive  anywhere ,  and  it  is  stable  if  Dp  is  negative 
everywhere. 

Adding  the  physical  requirement  that  w  must  be  con¬ 
tinuous  everywhere  to  (2.9),  we  see  that  if  p  and  Dp  are 
continuous  then  D(pDw)  must  also  be  continuous.  This  in 
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turn  implies  that  Dw  must  also  be  continuous.  Thus  the 
continuity  of  p  and  Dp  implies  the  continuity  of  w  and  Dw. 

If  p  or  Dp  is  discontinuous  (say  at  z  =  0) ,  then  a 
boundary  condition  can  be  derived  by  integrating  (2.9)  over 
an  infinitesimal  interval  including  2=0.  The  result  is 

Aq  (pDw)  =  -  (AoP)  w  (0)  (2.10) 

where  A0  \p  =  \p(0+)  - 

For  the  hydromagnetic  case  with  a  magnetic  field  in  the 
x-y  plane  equation  (2.9)  becomes  (Ref  2:55) 

D{  [py2  +  (&  •  S)2]  Dw 

=  k2  [py2  +  ~  (£  •  S)2  -  g  Dp]  w  (2.11) 

This  equation  is  more  useful  when  transformed  into  a  varia¬ 
tional  form.  This  is  done  by  multiplying  (2.11)  by  w  and 
integrating  both  sides  over  some  interval  of  interest 
(a,b) ,  giving 

/bw  D{[py2  +  -fi  (Ic  •  S)2]  Dw}  dz 

cl  "  41T  " 

*  k2  /b  w2[py2  +  TT  *  S)2  -  gDp] dz 

(2.12) 

The  left-hand  side  of  (2.12)  is  then  integrated  by  parts, 
assuming  the  boundary  condition 

w  Dw  [py2  +  (ic  •  $)2]  b  =  0  (2.13) 

a 
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This  condition  is  met  in  all  cases  considered  here  and 


(2.12)  then  becomes,  after  rearranging, 

2  {{Dp)w2  — [(Dw)2  +  k2w2]  }dz 

= 5 - 3— 42 -  (2-14) 

af  p{(Dw)2  +  k2w2}dz 

This  is  the  variational  form  which  will  be  used  here. 
Examination  of  (2.14)  shows  that  it  gives  a  real  value  for 
Y2  if  k  is  real  (which  it  is  here) .  Also,  since  the 
integral  of  the  second  numerator  term  in  (2.14)  is  positive 
definite,  we  see  that  the  growth  rate  decreases  when  a  mag¬ 
netic  field  is  introduced  with  a  component  parallel  to  the 
perturbation  wavevector  k.  This  stabilization  is  due  to 
the  supporting  magnetic  pressure  Pmag  which  helps  to  counter 
the  plasma  acceleration  g.  If  ic  1  £,  the  stabilizing  term 
vanishes. 

Rayleiqh-Taylor  Instability 

The  Rayleigh-Taylor  instability  is  a  macroscopic 
effect  often  encountered  in  hydrodynamics  and  hydromagnetics. 
It  occurs  when  a  more  dense  fluid  is  supported  by  a  less 
dense  fluid  against  a  "downward"  force.  A  classic  example 
is  a  layer  of  water  supported  by  a  layer  of  oil  against 
gravity.  A  small  periodic  perturbation  along  the  fluid 
interface  initially  grows  exponentially  in  amplitude.  This 
growth  can  be  explained  with  a  simple  energy  argument. 

Figure  3  shows  the  perturbation.  This  perturbation  involves 
the  transfer  of  fluid  from  positive  z  to  negative  z  regions. 
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This  is  an  energetically  favorable  transfer.  The  per¬ 
turbation  will  therefore  grow.  It  has  been  shown  experi¬ 
mentally  (2:51)  that  this  exponential  growth  changes  into 
a  linear  growth  with  respect  to  time  when  the  amplitude  of 
the  perturbation  reaches  about  2.5/k.  In  this  stage  spikes 
of  the  heavier  fluid  fall  and  bubbles  of  the  lighter  fluid 
rise  across  the  interface.  Paradoxically  the  initial  ex¬ 
ponential  growth  is  called  the  linear  Rayleigh-Taylor  in¬ 
stability  while  the  linear  growth  stage  is  called  the  non¬ 
linear  Rayleigh-Taylor  instability  (because  the  linear  MHD 
equations  no  longer  apply) . 

To  show  the  linear  instability  analytically,  consider 
the  Sturm-Liouville  equation  (2.9).  This  equation  has 
simple  analytical  solutions  for  two  different  density  dis¬ 
tributions.  Rayleigh  (15)  showed  that  for  an  exponential 
density 

P  =P0eez  (2.15) 

trapped  between  rigid  walls  at  z  =  0  and  L,  the  velocity 
perturbation  is  of  the  form 

w  =  A  (eml2  -  em2z)  (2.16) 

where  m^  and  m2  are  the  roots  of 

+  me  -  k^  (I-gB/y^)  =  0  (2.17) 

with  k,  g,  and  y  having  their  previous  meanings.  The 
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Rayleigh-Taylor  instability  mechanism. 
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instability  growth  rate  is  then  the  solution  of  the  rela¬ 
tion 


24  =  1  +  T  *2h2+*2  (2.18 

77  — vJY?— 

Since  the  right-hand  side  of  (2.18)  is  positive,  y^  takes 
on  the  sign  of  8.  The  fluid  is  therefore  stable  for  nega¬ 
tive  8  and  unstable  for  positive  8. 

Taylor  (18)  considered  the  stability  of  a  discon¬ 
tinuous  density  of  the  form 


for  constant  p.^  and 
w(z) 


p  =  PX  ,  Z  <  0 
p2  /  2  >  0 

P2*  Equation  (2.9) 

=  Ae^2  ,  z  <  0 
Ae-kz  ,  z  >  0 


(2.19) 

then  has  a  solution 

(2.20) 


Pl  <  p2  is  thus  an  unstable  configuration,  while  p-^  >  P2 
is  stable. 

This  section  has  dealt  so  far  only  with  the  hydro- 
dynamic  case.  There  is  a  completely  analogous  effect  in 
the  hydromagnetic  case  in  which  the  pressure  P  in  (2.1) 
is  replaced  by  a  sum  of  plasma  and  magnetic,  pressures 
Pmat  +  Pmag*  The  nature  of  the  instability  is  basically 
unchanged,  however. 
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A  Continuous  Density  Transition  (14) 

The  linear  growth  rate  estimates  currently  used  in  the 
SHIVA  project  are  due  to  LeLevier,  Lasher,  and  Bjorklund 
(14).  They  extended  Taylor's  treatment  of  a  discontinuously 
stratified  medium  by  assuming  a  continuous  density  transi¬ 
tion  of  the  form 


P  “  Pi  -  \  6Pe~kz  ,  z  >  0 


P2  +  6pe^z 


,  z  <  0 


(2.22) 


with  Sp  =  pi  -  P2  >  0.  Even  for  this  relatively  simple 
density  distribution  the  differential  equation  (2.9)  has 
no  obvious  analytical  solution.  LeLevier  et  al  instead 
begin  with  the  plasma  equation  of  motion 


+  DP  +pg=0  (2.23) 

and  the  assumption 

P1  =  -/wDpdt  (2.24) 


for  the  perturbed  density.  (2.23)  and  (2.24)  are  equiva¬ 
lent  to  the  hydrodynamic  equations  (2.1)  -  (2.6)  for 
B  =  0. 

A  velocity  perturbation  of  the  form 


w  =  Af (t)  cos  kx  exp  (-kz) ,  z  >  0 
Af (t)  cos  kx  exp  (+kz) ,  z  <  0 


(2.25) 


is  assumed,  presumably  because  (2.25)  is  the  simplest  solu¬ 
tion  of  the  incompressibility  condition 


V  *  $  =  0 


(2.26) 


in  rectangular  coordinates.  (2.23)  is  then  integrated 
over  the  intervals  (-°°,0)  and  (0,  =°)  using  the  velocity 
(2.25).  This  gives  two  expressions  for  the  pressure  at 
z  =  0  which  then  are  equated.  Quadratic  terms  in  f  and  f 
are  small  and  are  dropped.  The  resulting  equation  is 
differentiated  with  respect  to  time  to  give 

n  VK 

(P1  +  p2)f  "  9  fcfK  «Pf  =  0  (2.27) 

The  solution  of  this  is 


f(t)  =  exp  ( yt) 


(2.28) 


with 


2  =  gkK  6p 
Y  (k+K) (P1+P2) 


(2.29) 


An  interesting  characteristic  of  this  growth  rate  estimate 
is  that  y  remains  finite  for  short  wavelengths  (k  -*■  ®)  . 

It  also  agrees  well  with  numerical  simulations  of  SHIVA 
implosions  (17) . 

The  assumed  velocity  (2.25)  is  not,  however,  consis¬ 
tent  with  the  boundary  condition  (2.10)  or  the  differential 
equation  (2.9)  .  There  is  no  better  choice,  however,  that 
satisfies  the  incompressibility  condition  (2.26).  This 
problem  has  no  obvious  solution. 

It  should  be  mentioned  that  LeLevier  et  al  have  two 
sign  errors  in  their  report.  First,  the  gravitational 
term  in  (2.23)  has  the  wrong  sign  in  their  report.  Second, 
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the  provided  velocity  potential  $  which  leads  to  (2.25)  by 
the  relation  w  =  — D4>  requires  a  negative  sign  for  the 
negative  branch. 

Linear  Rayleigh-Ritz  Method  (6:193-206,  12:161-4) 

The  Rayleigh-Ritz  method  is  a  simple  variational 
method  used  to  approximate  the  eigenvalues  of  an  ordinary 
differential  equation.  Only  its  application  to  the  Sturm- 
Liouville  equation 

D(p  Dw)  -  rw  +  Xsw  =  0  (2.30) 

where  p,  r,  s,  and  w  are  all  functions  of  z  and  all  eigen¬ 
values  are  positive,  will  be  discussed  here. 

As  shown  in  Gelfand  and  Fomin  (6:198),  (2.30)  can  be 
expressed  in  the  variational  form 

/b  (p(Dw)2  +  rw2)dz 

X=  - - — -  (2.31) 

rD 

aJ  sw2  dz 

If  both  numerator  and  denominator  are  positive  definite 
quantities,  Rayleigh-Ritz  theory  states  that  the  function 
w  satisfying  appropriate  Boundary  conditions  and  minimizing 
X  in  (2.31)  is  an  eigenfunction  of  equation  (2.30).  This 
eigenfunction  corresponds  to  the  smallest  eigenvalue  of 
(2.30).  This  eigenvalue  is  then  given  by  (2.31)  with  the 
eigenfunction  substituted  for  w.  Normally,  however,  this 
method  is  used  only  when  the  eigenfunctions  cannot  be 
found  explicitly. 
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A  trial  function  with  adjustable  parameters  is  used  to 
approximate  the  eigenfunctions.  In  a  linear  Rayleigh-Ritz 
approximation,  the  trial  function  is  of  the  form 

q 

w  =  (z)  (2.32) 

where  the  <j>'s  are  orthonormal  functions  satisfying  the 
necessary  boundary  conditions.  This  trial  function  is  now 
inserted  into  equation  (2.31),  giving  an  estimate  of  the 
eigenvalue.  According  to  Rayleigh-Ritz  theory  (6) ,  this 
estimate  is  greater  than  the  true  eigenvalue  and  converges 
monotonically  to  the  eigenvalue  as  q  increases. 

The  problem  is  now  restated  as  that  of  minimizing 

J  (aj^,  ...»  ctq)  =  f  (p(Dw)2+rw2)dz  (2.33) 

subject  to  the  constraint 

K  (aj_,  ...,  a_)  =  j-**  sw^dz  =  1  (2.34) 

H  a 

This  problem  is  identical  to  the  original  problem  and  in  no 
way  limits  the  class  of  allowable  solutions  (6) .  The  new 
problem  is  solved  using  Lagrangian  multipliers.  Minimiz¬ 
ing  J  subject  to  K  =  1  therefore  requires  that 

-r— —  (J  -  oK)  =  0  ,  i=l ,  ...,  q  (2.35) 
d  a 

where  a  is  an  undetermined  Lagrangian  multiplier.  (2.35) 
is  a  system  of  q  simultaneous  linear  equations.  If  Ji  j 
and  K^j  are  defined  such  that 


19 


q  q 

J  “  ill  j =i  ai°j  Jij 

(2.36) 

and 

q  q 

K  “  i=l  i=l  ai“j  Kij 

(2.37) 

then  the 

system  will  have  a  nontrivial 

solution  if  and  only 

if 

det  (D)  -  0 

(2.38) 

where 

. 

Di j  “  Jij  “  a  Ki j 

(2.39) 

o  is  then  the  best  estimate  for  the  smallest  eigenvalue 
A^  of  (2.30)  possible  with  the  trial  function  (2.32)  of 
order  q.  As  the  number  of  terms  q  in  the  trial  function 
increases,  the  estimates  decrease  monotonically  towards 
A^  as  q  goes  to  infinity. 

The  Rayleigh-Ritz  method  thus  provides  an  upper  bound 
on  the  first  eigenvalue  Aj  of  the  Sturm-Liouville  problem. 
The  number  of  terms  in  the  trial  function  can  be  increased 
until  the  accuracy  of  the  approximation  is  sufficient. 
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III.  Applications 

Application  of  Rayleigh-Ritz  Theory  to  the  Hydromagnetic 
Problem 

The  Rayleigh-Ritz  theory  of  Chapter  II  will  here  be 
applied  to  the  variational  equation  (2.14)  to  estimate 
the  growth  rate  y  for  an  accelerating  plasma  sheath.  The 
sheath  density  will  be  assumed  to  peak  in  the  center  of 
the  sheath  with  Dp  >  0  over  0  <  z  <  zQ  and  Dp  <  0  for 
z  >  z0  (see  Figure  4) . 

>  The  maximum  growth  rate  might  be  expected  with  a  trial 
velocity  function  w  which  is  identically  equal  to  zero  for 
z  >  zQ.  This  assumption  implies  that  there  is  no  net  plas¬ 
ma  flow  from  the  stable  (Dp  <  0)  into  the  unstable  (Dp  >  0) 
region.  As  shown  in  Chapter  II,  Dw  must  be  continuous 
everywhere  if  p  is  well  behaved.  Both  w(z0)  and  Dw(z0) 
must  therefore  equal  zero. 


Figure  4.  Example  of  the  class  of  density 
profiles  treated  here. 
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An  orthonormal  set  of  functions  is  needed  for  the 


Rayleigh-Ritz  procedure.  The  functions  must  be  defined 
over  the  interval  (0,  zQ)  and  meet  the  boundary  conditions 
w(zQ)  =  Dw  (z0)  =  0.  Functions  satisfying  these  conditions 
are 


where 


cm  ~  |)  and 


Sm 


and 


cosh  (Xmu)  _  cos  ( Xmu) 
cosh  (Xm/^)  cos 


Sm(u) 


sinh  (ymu)  _  sin  (umu) 
sinh  (um/^)  sin  (ym/^) 


(3.1) 

(3.2) 


Cm  and  Sm  are  discussed  in  Appendix  B.  They  satisfy  the 
boundary  conditions 

Cm  (±  \)  =  (±  J)  =  (±  |)  *  Sm  (±  |)  *  0  (3.3) 

\  satisfactory  trial  velocity  function  is 

=  m=l^amCm  (z£  l)  +  emSm  2)  ^  (3.4) 

(3.4)  also  satisfies  w(0)  =  Dw(0)  =  0.  These  boundary 
conditions  are  equivalent  to  placing  a  rigid  wall  at 
z  =  0.  For  lack  of  a  better  trial  function,  (3.4)  will 
be  used.  Squaring  (3.4), 

”2  '  mil  nil  '“mCm  <if  -  |>  +  -  §) ) 

•  («nCn  'if  -  l’  +  8nS"<it  '  l"  (3-5) 
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Differentiating  (3.4)  and  then  squaring  gives 


M2  -  TJ  mil  nil  -  |) ) 


(°ncn  (Zq  -  2^  +  6nsn<Zo  ”  2^ 


(3.6) 


According  to  the  linear  Rayleigh-Ritz  theory  of  Chapter  II, 
the  required  equations  are 


and 


3  3m 


3  3m 


(J  -  XK)  =  0  m=l , 


(J  -  XK)  =  0  m=l ,  ...,  q 


(3.7) 


(3.8) 


where  K  and  J  are  the  numerator  and  denominator,  respec¬ 
tively,  of  (2.14).  X  is  the  Lagrangian  multiplier. 

Inserting  (3.5)  and  (3.6)  into  (2.14),  we  obtain 

Z° 

J  =  /  p((Dw)2  +  k2w2)dz 


z0 

1  p{ 


1 

Zq' 


mil  nIl(amV 


Cm(~  - 


2*  +  6msm(2Q  “  2^ 


•  *°ncn  'if  -  ?»  +  6nsn  'if  ‘  » 

+  k2  mil  nil '“-"Cm 'if  -  ♦  "mSm'if  -  §> ) 

•  '°ncn'if  -  +  Snsn(ji  -  ill)  dz  0.9) 

Making  the  substitution 
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using  the  notation 


(♦llf  1 4»2)  =  /  $l(u)  f  (u)  (j>2  (u) 
and  differentiating  with  respect  to  am, 

^  =  ~  nli  U°n(Cm|p|cn)  +  en  <Cni|p|S£)) 

+  k2z02  (an (Cm| p | Cn)  +  en(Cm|p|Sn)) }  (3.10) 


Noting  the  symmetry  between  C  and  S  in  (3.9),  we  can  simply 
interchange  a  with  @  and  C  with  S  in  (3.10)  to  obtain 


3  J 
3  em 


nil  f  ^“n  (cr 


&n  ^  sm  I 


+  k2Zo2  (otn  (Cj  p  1%)  +  en(SnJp|Sn))  }  (3.11) 


The  function  K  will  be  subdivided  into  K  =  -  K2r 

with 

zo 

Ki  =  /  (Dp)w2  dz  (3.12) 

o 

and 

K2  =  Sgsifi-  /Z°  B2  {(Dw)2  +  k2w2 }  dz  (3.13) 

4ug  Q 

Performing  exactly  the  same  operations  on  and  K2  as 
were  performed  on  J  gives 


fU  *  2zo  nli  (an  (Cm|  Dp  | Cn)  +  Bn  (Cm|  Dp  |  Sn)  )  (3.14) 

q 

ffi  -  2zo  nIx  (an  (Cn !  Dp  |  Sm)  +  6n(Stn|Dp|Sn))  (3.15) 
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3K2  _  2  cos2  9 
8am  4ugzo 


nSiUantCmlBZlCfl)  +  fin  (C£|  b2  |  S„) ) 


+  k2Zo2  (an  (Cj^I  B2  |  Cn)  +  en(Cm|B2|sn)  )  }  (3.16) 


and 


3K2  _  2  cos2  e 
3  0m  4irgz0 


{ (an(Cn|B2|s^)  +  Sn(S^|B2 | S^) ) 


+  k2z02  (anlCnlBZlSm)  +  Bn (Sm|  B2 | Sn) ) }  (3.17) 

To  simplify  the  final  result,  the  relations  (3.7)  and  (3.8) 
can  be  redefined  as 

'aTT^  ^  *  n-1  ^amn  “n  +  ^mn  ®n)  "  ®  (3.18) 

and 

-  (J  -  AX)  =  (cmn  an  +  dmn  0n)  *  0  (3.19) 

Once  p  and  B  have  been  inserted  in  (3.10),  (3.11),  and 
(3.14)  -  (3.17)  the  elements  a^,  bum,  c^,  and  dmn  can  be 
determined.  With  these  new  variables  the  equation  to  be 
satisfied  is  now  in  the  block  determinant  form 


A 


C 


(3.20) 


with  A,  B,  C,  and  D  being  q  by  q  matrices. 

Once  the  matrix  elements  of  (3.20)  have  been  de¬ 
termined  for  a  particular  problem  the  values  of  Az0  (a 
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dimensionless  quantity)  which  solve  (3.20)  can  be  found 
as  a  function  of  k.  Applying  the  Rayleigh-Ritz  theory 
from  Chapter  II,  we  see  that  the  growth  rate  estimate  for 
(2.4)  is 


or 


_li  a 
W2 


1 

opt 


y2  =  jc?_8Q2  _  _2 

^OptZo  zo 


(3.21) 

(3.22) 


where  the  optimal  root  A0pt  is  the  smallest  positive  root 
of  (3.20)  or  the  largest  negative  one  if  there  is  no 
positive  root.  A0pt  therefore  corresponds  to  the  most 
unstable  mode  of  the  instability. 


Approximate  Forms  for  B  and  £ 

In  a  SHIVA  implosion,  the  magnetic  field  is  initially 
located  totally  outside  the  plasma  sheath.  The  field  then 
diffuses  into  the  sheath  according  to  the  diffusion  equa¬ 
tion 


3B  _  n  d2B 
3t  4 tt  TzZ 


(3.23) 


in  a  planar  approximation,  where  n  is  the  electrical  re¬ 
sistivity.  Hussey  and  Roderick  (10)  have  shown  that  the 
external  field  B0  increases  linearly  with  time  through 
much  of  the  implosion  time.  This  leads  to  a  magnetic 
field  profile 


B (z, t)  =  4  B0(t)  i2  erfc  (^)  (3.24) 
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where 


q  = 


(3.25) 


and 

i2  erfc  (^)  (3.26) 

is  the  second  repeated  integral  of  the  complementary  error 
function. 

The  greatest  contribution  to  the  instability  growth 
rate  comes  from  the  outer  edge  of  the  plasma  sheath  where 
the  plasma  density  rises  quickly.  For  small  z  (3.26)  is 
well  approximated  by  (10) 


where 


B  =  B0  (1  -  1  J)  (3.27) 

it  2  q 

q  =  4.51  Xth  (3.28) 


for  a  sheath  thickness  X^h*  Sheath  thickness  is  here 
defined  as  sheath  mass  per  unit  area  divided  by  maximum 
sheath  density. 

The  linear  magnetic  field  (3.27)  equals  zero  when 

1 

_  2 

z  =  =  2  Xth  (3.29) 

This  small  z  approximation  will  be  extended  throughout  the 
whole  plasma  sheath  under  the  assumption  that  the  small-z 
region  determines  the  instability  of  the  sheath.  The  mag¬ 
netic  profile  used  here  is  therefore 

B (z)  =  Bo  (1  -  f)  (3.30) 


27 


where  L  =  2Xtjj*  This  profile  is  illustrated  in  Figure  5. 

Hussey,  Roderick,  and  Kloc  (11:1454-5)  have  used  the 
equilibrium  equation  of  motion 

pg  =  “  3Y  (pmat  +  pmag)  (3.31) 

to  derive  an  equilibrium  density  distribution  which  cor¬ 
responds  to  the  linear  magnetic  field  profile  (3.30).  This 
density  distribution  is  (see  Figure  6) 


p 

% 


2p  (1  +  6  -  |  -  (1+8)  exp  (-  ~) )  ,  0  <  z  <L 

2p  (6  -  (1+8)  exp  (-  ^) )  exp  (^^)  ,  z  >L 

p  pL 


(3.32) 


where  p  is  an  average  sheath  density  and 

_  pmat 
pmag 


(3.33) 


is  defined  in  terms  of  the  maximum  magnetic  pressure  and 
the  average  plasma  pressure.  Differentiating  p  with  respect 
to  z  shows  that  the  peak  density  occurs  at 


z  =  zQ  =  8L  In(^)  (3.34) 

and  is 

Pmax  =  2p  (1-8  Ln(±±S))  (3.35) 

For  small  values  of  8,  very  little  plasma  is  located  in  the 
region  z  >  L .  The  plasma  is  therefore  concentrqted  in  a 
region  twice  the  defined  sheath  thickness. 
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Figure  5 
magnetic 


max 


Growth  Rates  for  a  Plasma  Sheath 

Growth  rates  for  the  Rayleigh-Taylor  instability  in 
an  accelerating  plasma  sheath  are  addressed  in  this 
section.  This  development  is  the  major  result  of  this 
work. 

The  density  distribution 

p  =  2p  {(1+3)  -f  -  (1+b)  exp  (-  ~)  }  (3.35a) 

L  PL 

for  0  <  z  <  L  will  be  used  here,  along  with  the  magnetic 
field 

B (z)  =  B0  (1  -  f)  (3.36a) 

Ii 

(3.35a)  and  (3.35b)  need  only  be  inserted  into  the  equations 
(3.10)  -  (3.17).  The  determinant  equation  (3.20)  will  then 
be  solved,  giving  the  desired  growth  rate  estimate. 

If  the  density  (3.35a)  is  inserted  into  (3.10)  we  get 


3  J 
3otm 


4p 

zo 


(1  + 


-  a> 


q 

z 

n=l 


Uan(Cm|C')  +  Bn  (CrfilSn)) 


+  k2zQ2  (an(Cin|Cn)  +  &n^m|Sn))} 

~  l  { (an  (Cm | u | C^)  +  3n  (Cm| u [  Sn)  ) 


+  k^zQ^  ( <*n  (C^  |  u  |  Cn)  +  6n  (Cm  |  u  |  Sn)  )  } 

_  4~P  il±i>  exD  / _ Zfl) 

zQ  exp  '  2BL' 

*  Jl  i(«n(Cm'|exp(-  |CH)  +  8n<Cm|exp(-  |^)|S„)) 


+  k2Zo2(an(Cm|exp(-  |^)|Cn)  +  6n(Cm!exp(  -  |^)|Sn))} 


BL 


(3.36) 
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Orthonormality  yields 


(Cm  I  Cn)  —  and  (Cp^|sn)  —  0 


while  parity  considerations  mean  that 


(Cm I Sn)  =  (Cm | u | C£)  =  (Cm|u|Cn)  =  0 


Inserting  thsse  values  into  (3.36)  gives 

=  if  {1  ;  0  "  Il)  E.  an((Cm|Cn)  +  k2z02  6mn) 


-  4f  ?  en((Cm|uiS^)  +  k2z02(Cm|u|Sn)) 


L  n=l 


4p  (1+B) 


exP  (“  Jn.] 


•  J^ant^iexpt-  |Crf)  +  k2z02  (Cm|  exp  (-  ^)Cn)) 

+  6n(  (Cm'iexp(-  -^)|Sri)  +  k2zQ2  (Cm|  exp  (-  ^£)  Sn) )  } 


(3.36a) 


Noting  the  symmetry  between  C  and  S  in  (3.9),  we  can  simply 
interchange  a  with  6  and  C  with  S  in  (3.36a)  to  obtain 


^  (1  +  8  “  2*t.)  E.  ^n((Sm|Sn)  +  k2z02  6mn) 


2l'  nil 


-  i£  E 


nIl  an((Cn|u|Sm)  +  k2z02 (Cn | u | Cm) ) 


- ^  «p  <- 
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1 

i 


q 

*nIl{en<<Sm|exp(-  I  S£)  +  k2zQ2 (sm| exp (-  ^g)  |  Sn) ) 

+  an((Cn|exp(-  ^|)  |  S„)  +  k2Zo2 (cn ) exp (-  ^g)|Sn))}  (3.37) 

Inserting 

D»  *  2»  <-  l  +  TT  exP«-  sfl)  <3-38' 

into  (3.14)  and  (3.15)  yields 


3Ki  4pzo  ,  ,  1+6 

+ 


b  exp  '■  He' 


and 


*  n£i<«n<cm|exp(-  ^g)  |Cn)  +  Bn(Cm|exp(-  ^oH)  j  Sn) )  ) 

(3.39) 


3Kl  _  4pz0  {-6m  +  exp(-  |^f) 
L 


3  Bm 


2  6  L 


n£l  ^Bn (Sm| exp  (-  ~ 1 ^"H^n) 


+  an(Cnjexp(-  ^H)  |  Sm)  )  } 


(3.40) 


Inserting  the  B  field  (3.30)  into  (3.16)  and  (3. 17)  yields 


3K2  _  2Bq2ZqCOs2  e 
3  am  4irgL^ 


{an((C^|u2|Ch) 


+  k2Zo2(Cm|u2|cn) 


+  -  |)2((Cm|C')  +  6^  k2Zo2)) 

“  Bn  (=-  “  1)  ( (Cjn|  u|  SfJ)  +  k2Z{;)2  (C^l  u  |  Sn) )} 


(3.41) 


30m  °  4°gl^  (Bn(  (smlu2|Sn)  +  k2Zo2  (Sm|  u2  |  Sn) 

+  "  l)2((SmlSn)  +  k2zo2  W> 

-  «n<f"  D((Cn|u|Sm)  +  k2z02  (Cn  |  u  |  Sm) )  }  (3.42) 

The  expressions  (3.36)  -  (3.42)  can  now  be  substituted 
into  (3.18)  and  (3.19)  to  give  values  for  the  matrix  ele¬ 
ments  a'  through  d',  where 

L 

amn  =  p  amn  (3.43) 


and  so  on  for  b,  c,  and  d.  The  results  are  then 


aicin  **  x(l+8”  2^  ^cmlcn^ 

-  li^±§1  exp  (-  j|>  <<C*|®xp<-  iH)  |Cfi> 

+  k2z02  (Cm  |  exp  (-  ^)|Cn)) 

-  4xXL(^)  exp  (-  f|)  (Cm|exp(-  ^)|Cn) 

+  2TTXL(k2z02(Cm|u2|Cn)  +  (Cm|u2 |Cn) 

+  (7  -  7) 2  (cm|Cn)) 

+  6mn(4  U+8  -  J)  k2z02  +  4tXl+  2t  TXL  (£  -  J)2k2z02  ) 

P  *4) 


33 
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bmn  ~  cnm  ~  “4  (  (Cml  u  |  Sfj)  +  k2zo2  (Cm|u|Sn)j 


-  4  exp  (-  jg)  {(C^iexp(-  IH)  |s>)+k2Zo2(Cm|exp(-  I  Sn)  } 

-  4rXL(i±i)  exp  (-  (Cm|exp  (-  I|)  |  Sn) 

-  2tTXL(^  -  1)  {k2Zo2  (Cm|u|Sn)  +  (C^|u|s^)}  (3.45) 

and 

^mn  =  ~(l+3  -  2")  (Sm  |  S^) 


-  -y-]  exp  (-  ^ )  { ( Sm |  exp  (-  |Sn) 


+  k2Zo2 (sm|exp(-  ^|)|Sn)} 

-  4xXL(ij^)  exp  (-  j|)  (Sm|exp  (-  ^)|Sn) 


+  -  f) 

+  4tA  L  +  2TTJLk2z02(i  -  j)2} 

(3.46) 

with 

(3.47) 

and 

T  = 

Bo2cos2  e 

4TrgLlT 

(3.48) 

If  the  substitution 

_  _  Bq2 

^  -  8  ir  p  L 

(3.49) 

is  made,  then  (3.48) 

becomes 

T 

=  2  cos2e 

(3.50) 

All  matrix  elements  of  the  equation 


A' 


(3.51) 


are  now  known.  From  the  first  section  of  this  chapter, 
the  eigenvalues  X  of  this  equation  allow  an  estimate  of 
the  growth  rate 


or 


2  -  k2L2  2 
XoptL  L 


(3.52) 


v2L  _  k2Zp2 
g  Xopt  Lt2 


(3.53) 


A  FORTRAN  program  was  written  to  solve  (3.51)  for  X  as  a 
function  of  k2z02,  8,  and  8.  The  results  of  the  program 
are  discussed  in  Chapter  IV.  All  the  unevaluated  integrals 
in  (3.44)  -  (3.46)  are  evaluated  in  Appendix  C. 
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IV.  Results 


The  Rayleigh-Ritz  approximation  to  the  hydromagnetic 
instability  growth  rates  of  an  accelerating  plasma  sheath 
was  developed  in  Chapter  III.  A  computer  program  was 
written  to  implement  the  approximation.  The  results  of 
this  approximation  are  presented  here.  The  stabilizing 
effects  of  a  magnetic  field  diffused  into  the  sheath  are 
also  discussed. 

Parameters  needed  for  the  approximation  are  provided 

*w 

by  a  one-dimensional  simulation  to  the  SHIVA  implosion 
using  the  MAGPIE  code.  The  simulation  data  at  selected 
times  is  listed  in  Table  1.  This  raw  data  is  converted 
into  necessary  parameters  in  Table  2.  The  necessary  rela¬ 
tions  to  generate  these  parameters  are  given  in  Appendix  A. 

Convergence  of  the  Approximation 

As  the  number  of  terms  in  the  trial  velocity  w  in¬ 
creases,  the  growth  rate  approximation  provided  by  the 
procedure  of  Chapter  III  should  converge  to  the  actual 
growth  rate.  This  is  illustrated  in  Figure  7  for  several 
cases  with  perturbation  wavevector  ic  .1  2.  For  these  cases 
the  growth  rate  estimates  increase  by  25  to  30  percent  as 
the  trial  function  increases  from  2  to  4  terms.  The  growth 
rates  increase  by  just  1  to  2  percent  as  the  trial  function 
increases  from  18  to  20  terms.  The  convergence  seems  to 
be  rather  slow. 
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Instabilities  with  ic  l  5 

The  limiting  case  of  a  perturbation  £  J_  S  will  be  dis¬ 
cussed  here.  This  planar  case  is  equivalent  to  a  cylindrical 
instability  with  only  z-axis  perturbations.  In  cylindrical 
coordinates  it  is  just  the  sausage  instability. 

Figure  8  shows  the  resulting  dispersion  relation  for 
typical  plasma  parameters.  We  see  that  the  sheath  is  un¬ 
stable  for  all  wavelengths  of  perturbation,  approaching 
zero  for  long  wavelengths  and  an  upper  limit  for  the  growth 
rate  as  \  +  0  (k  -*•  ®)  .  This  upper  limit  is  not  evident  in 
the  figure,  but  it  was  verified  for  extremely  small 
wavelengths  (k  zQ  =ic£) .  This  general  behavior  agrees  with 
that  of  Tsai,  Liskow,  and  Wilcox  for  basically  the  same 
problem  (Figure  9) . 

The  dimensionless  growth  rates  of  Figure  8  are  trans¬ 
formed  into  actual  growth  rates  in  Figure  10  for  several 
times  during  the  SHIVA  implosion.  The  data  from  Table  2 
has  been  used,  with  the  8  rounded  to  the  nearest  multiple 
of  .05.  The  solid  lines  in  Figure  10  represent  growth  rate 
estimates  made  using  LeLevier's  result  for  an  exponentially 
transitioning  density.  It  is  assumed  here  that  LeLevier's 
density  scale  length  1/K  =  Zo/4ir.  The  agreement  between 
the  two  methods  is  then  excellent  (within  5  percent) . 

Oblique  Perturbations 

Perturbations  with  wavevectors  not  perpendicular  to  B 
were  considered.  Figure  11  shows  growth  rates  for  a 
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constant  perturbation  wavelength  as  the  angle  6  between 
it  and  5  decreases  from  90°.  As  the  wavevector  rotates, 
the  destabilizing  perpendicular  component  decreases  while 
the  stabilizing  parallel  component  increases  (parallel  and 
perpendicular  will  hereafter  be  with  respect  to  §) .  There 
therefore  exists  an  angle  6marg  at  which  marginal  stability 
is  achieved.  This  ©marg  is  an  overestimate  which  also 
converges  to  the  true  value  as  the  length  of  the  trial 
function  increases. 

s  The  variation  of  6marg  with  perturbation  wavelength 
and  the  number  of  trial  function  terms  is  illustrated  in 
Figure  12.  We  see  that  0marg  increases  moderately  as 
perturbation  wavelength  decreases.  Even  though  the  un¬ 
stabilized  perpendicular  instability  is  greater  for  larger 
wavenumbers,  the  stabilizing  influence  of  the  parallel 
instability  grows  even  faster  with  3'.  Stability  therefore 
occurs  at  smaller  angles  as  k  increases. 

In  Figure  13  the  perpendicular  component  of  the  per¬ 
turbation  wavevector  kx  remains  constant  while  the  parallel 
component  ky  is  increased.  In  this  case  the  destabilizing 
effect  of  kx  is  held  constant  while  the  stabilizing  influ¬ 
ence  of  ky  is  gradually  increased  until  a  value  (ky)marg 
is  reached  at  which  the  two  effects  cancel  out  to  yield 
marginal  stability.  Figure  13  shows  that  long-wavelength 
parallel  components  can  be  used  to  stabilize  much  shorter- 
wavelength  perpendicular  components  \x,  with  the  growth 
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rate  being  linear  in  ky.  This  stabilization  causes  the 
growth  rates  to  decrease  when  a  long-wave length  component 
is  introduced.  This  fact  is  especially  important  when  the 
problem  is  related  to  cylindrical  geometry.  The  stabilizing 
9-perturbation  wavelengths  of  a  cylindrical  sheath  are  on 
the  order  of  the  cylinder  radius  for  low-order  modes  while 
the  destabilizing  z-perturbation  wavelengths  are  on  the 
order  of  the  sheath  thickness.  This  suggests  that  long- 
wavelength  9  perturbations  might  be  able  to  largely 
Stabilize  the  z  instabilities. 

These  results,  however,  disagree  with  previous  work 
done  on  this  problem.  Tsai,  Liskow,  and  Wilcox  (19)  have 
recently  reported  the  dispersion  relation  for  parallel 
perturbations  shown  in  Figure  14.  The  sheath  becomes  un¬ 
stable  for  long  wavelength  perturbations,  with  a  mode  of 
maximum  instability  existing  approximately  midway  between 
the  two  points  of  marginal  stability.  These  unstable  modes 
were  not  observed  here,  however. 

Experimental  evidence  also  supports  the  belief  that 
unstable  modes  exist  for  long  wavelength  perturbations 
parallel  to  §.  SHIVA  implosions  have  been  distorted  most 
seriously  by  azimuthal  perturbations  in  a  pentagonal  shape 
(ke  =  5) .  These  perturbations  have  wavelengths  on  the 
order  of  the  cylinder  radius. 

These  unstable  modes  might  not  be  observable  here 
because  the  number  of  terms  in  the  trial  function  needs 


to  be  increased.  As  the  length  of  our  trial  velocity 
function  increases,  the  growth  rate  estimate  increases. 

The  expected  parallel  instability  growth  rates  may  be  so 
small  that  short  trial  functions  will  indeed  lead  to 
negative,  stable  growth  rates.  Longer  trial  functions  would 
then  be  needed  to  give  positive  growth  rate  estimates. 


#  of  Trial  Function  Terms  (2q) 


Figure  7.  Dimensionless  growth  rate 
estimates  as  a  function  of  trial  func¬ 
tion  length. 


Figure  o.  Dimensionless  growth  rate 
estimates  as  a  function  of  perturbation 


wavenumber  for  0=90°. 


Figure  9.  Dispersion  relation 
of  Tsai  et  al  for  ©=90°.  (19) 
The  dimensionless  growth  rate 
is  plotted  as  a  function  of  the 
dimensionless  perturbation 
wavelength. 
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-  My  estimates  with  20  term  j 

trial  function  — ; - j— 

+  LeLevier  et  al  estimates  j 

0  20  40  60  80  100 


Figure  10.  Growth  rate  estimates  as  a  function 
of  perturbation  wavenumber  for  true  SHIVA 
parameters  and  ©=90°  compared  to  LeLevier  et 
al's  growth  rate  (2,29)  »  which  assumes  an 
exponential  density  transition. 


20  Term  Trial  Function  Used 
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Figure  13.  Dimensionless  growth  rate  versus  the 
wavenumber  component  parallel  to  the  magnetic  field, 
k^  is  the  component  perpindicular  to  the  magnetic 
field  while  k  is  the  component  parallel  to  the 
magnetic  field. 


Figure  14.  Dispersion  relation 
of  Tsai  et  al  for  0=0°,  (19) 

The  dimensionless  growth  rate 
is  plotted  as  a  function  of  the 
dimensionless  perturbation 
wavelength. 
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V.  Conclusion 


A  method  has  been  developed  to  approximate  the  Rayleigh 
Taylor  growth  rate  of  an  accelerating  plasma  slab.  It  has 
been  shown  to  agree  with  the  analytical  results  of  LeLevier 
et  al  (14)  for  the  limiting  case  ic  ]_  3.  This  agreement 
suggests  that  my  choice  of  a  trial  velocity  function 
identically  equal  to  zero  in  the  stable  region  is  valid. 

The  behavior  of  the  other  limiting  case,  ic  parallel  to  3, 
however,  is  suspect.  No  unstable  modes  were  found  for 
long  wavelengths.  This  fact  may  be  due  to  insufficiently 
long  trial  velocity  functions.  If  this  is  not  the  case, 
then  there  may  be  a  more  serious  problem  in  this  work 
which  brings  the  presented  results  (at  least  for  ic  not 
perpendicular  to  3)  into  question.  As  mentioned  in  the 
Introduction,  the  perturbations  parallel  to  3  are  much  more 
interesting  than  those  perpendicular  to  3. 

The  developed  program  can  also  be  used  to  estimate 
the  angle  between  ic  and  3  at  which  a  particular  perturba¬ 
tion  wavelength  becomes  marginally  stable.  For  a  given 
wavelength  the  instability  has  been  shown  to  stabilize 
♦rapidly  as  0  decreases  from  90°. 

The  method  used  here  is  general  and  can  be  used  for 
any  density  profile  with  only  one  contiguous  unstable 
region.  The  only  change  necessary  is  to  replace  p  and  B  in 
the  unevaluated  integrals  in  the  first  section  of  Chapter 
III.  Densities  with  several  unstable  regions  simply 
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require  a  trial  velocity  function  which  is  piecewise 
defined. 

This  work  can  be  continued  in  several  ways.  First, 
estimates  could  be  made  for  a  greater  variety  of  parameter 
values.  The  present  results  represent  a  fairly  narrow 
parameter  range.  In  particular,  longer  trial  velocity 
functions  need  to  be  used  with  long-wave length  perturbations 
parallel  to  B.  The  expected  instabilities  might  then  be 
observed.  Second,  the  method  used  to  find  the  solutions 
of  the  determinant  equation  (3.20)  —  essentially  a  binary 
search  method  —  might  be  improved  to  minimize  the  number 
of  times  the  determinant  must  be  evaluated.  The  current 
method  uses  a  prohibitive  amount  of  computer  time  for  long 
trial  functions  (and  therefore  large  determinants).  Third, 
more  accurate  density  and  magnetic  field  profiles  could 
possibly  be  found.  Hussey  and  Roderick  (10)  have  presented 
more  accurate  profiles  which  are  probsbly  too  difficult  to 
work  with  here.  Better  approximations  to  these  functions 
than  the  ones  used  here  might  be  possible,  however. 
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Appendix  A:  Cylindrical  Sheath  Implosion  Relations 


The  basic  variables  for  the  implosion  of  a  cylindrical 
plasma  sheath  can  be  related  by  a  series  of  simple  relations 
These  relations  are  useful  enough  to  be  derived  here. 

The  cylindrical  configuration  is  shown  in  Figure  15. 

The  first  equation  relates  the  magnetic  field  to  the  current 
passing  through  the  sheath.  Outside  the  sheath  the  field 
is  simply  that  of  a  thin  wire  on  axis: 


b0 


(G) 


=  .21  (A) 
r(cm) 


(A.  1) 


where  r  is  the  radial  coordinate. 

Equation  (A.l)  can  be  used  to  find  the  plasma 
parameter 


_  Pmat 
Pmag 


(A.  2) 


where  Pmat  is  the  average  material  pressure  and  Pmag  the 
maximum  magnetic  pressure  in  the  sheath.  We  may  take 


pmat  “  2  Ppeak 


(A. 3) 


where  Ppeak  is  the  peak  plasma  pressure  (which  is  the 
available  variable  from  computer  simulations)  and 


Pm  (Mbar) 


Be2  (rQ)  _  Bq2  (~MG) 

8ff  8  TT 


(A. 4) 


for  a  cylinder  radius  rQ 


The  next  variable  needed  is  the  sheath  acceleration  g„ 
It  is  given  by  the  relation 


g  = 


Bol 

8ito 


(A. 5) 


where  a  is  the  areal  mass  density  (7:1057).  If  a^n  is  the 
areal  mass  density  of  the  initial  cylindrical  foil,  then 
a  is  governed  by  a  simple  cylindrical  convergence  relation: 

ainrin 


a  = 


ro 


(A. 6) 


where  r^n  is  the  initial  foil  radius.  Inserting  (A. 6) 
into  (A. 5)  then  gives 


9  - 


Bo2rin 

®woinro 


(A.  7) 


The  final  relationship  will  be  for  the  sheath  thick¬ 
ness  Atft.  Hussey  and  Roderick  (9:1385)  have  developed 
the  approximate  relationship 


,  _  1  ,nt*  i 

Xth  “  4.5i  <  tt  2 


(A. 8) 


where  n  is  the  electrical  resistivity  and  t  is  the  time 
measured  from  the  beginning  of  the  implosion.  For  lack  of 
better  information  the  resistivity  will  be  assumed 
constant  here. 
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Figure  15.  Geometry  of 
cylindrical  implosion. 


•  ■  Appendix  B:  A  Set  of  Orthonormal  Functions 

A  set  of  orthonormal  functions  satisfying  the  boundary 
conditions 


♦  (a)  =  <T(a)  =  <t>  (b)  ■♦'(b)  =  0  (B.l) 

for  a  <  b  is  needed  in  Chapter  III.  The  functions 


C™<u>  =  co £ 

(Xmu) 

cos 

( *mu) 

(B.  2) 

(W2> 

cos 

and 

V 

s™<u>  -  flHiT 

(ymu) 

"W?) 

sin 

sin 

(wmu) 

Tum7^) 

(B.  3) 

were  discussed  by  Chandrasekhar  (3:634-643)  and  satisfy 
the  equation 

y(iv)  -  tt4y  (B.4) 

with  the  boundary  conditions 

y(±|>  =  y'(*§)  =  o  (b. 5) 

The.Cs  and  S's  are  orthonormal  over  the  interval  (-  ^ ,  ^-) 
with  eigenvalues  given  by  the  solutions  of 

tanh  j  +  tan  ^  -  0  (B.6) 

and 

coth  -  cot  Jjr  =  0  (B.7) 

The  functions  Cm  and  Sm  for  m  =  1,2,3,  and  4  (and 
their  first  three  derivatives)  have  been  tabulated  from 
u  =  0  to  0.5  by  Harris  and  Reid  (16).  Basic  indefinite 
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¥ 


1 


integrals  involving  Cjjj  and  Sm  have  been  published  by 
Reid  and  Harris  (16) .  These  results  allow  integrals  of 
the  form 


($1 1  vin  |  4>2 ) 


*5 

/  4>  1  (u)  un  4>2  (u)  du 


(B.8) 


to  be  calculated  simply,  with  and  <f>2  representing  C  or 
S  and  n  integer.  (4>i  |  un  j  412)  is  integrated  by  parts  q 
times  to  leave  only  integrals  of  the  form  |<J>2^S^)» 

The  indefinite  forms  of  these  integrals  are  evaluated  in 
Reid  and  Harris  for  r  =  0,1,2,  and  3.  Combining  these 
indefinite  integrals  with  the  boundary  conditions  (B.5) 
gives  a  value  for  the  integral  (B.8).  The  explicit  forms 
of  Cm  and  Sm  have  thus  been  avoided. 


59 


Appendix  C:  Evaluation  of  Needed  Integrals 


All  integrals  of  the  form  ( !  f  I  <^2 )  used  in  Chapter 
are  evaluated  here.  These  integrals  fall  into  two  cate¬ 
gories:  (1)  those  with  f (u)  =  1,  u,  or  u2  and  (2)  those 

with  f(u)  =  exp  (-au) .  Using  the  C  and  S  functions  of 
Appendix  B  as  <{>]_  and  $2*  the  first  category  of  integrals 
can  be  handled  fairly  easily  using  indefinite  integrals 
solved  in  an  article  by  Reid  and  Harris  (15) .  The  second 
category  of  exponential  integrals,  on  the  other  hand,  must 

s. 

be  integrated  directly  using  the  explicit  forms  (B.2)  and 
(B.3)  of  C  and  S. 

As  an  example  of  the  method  used  for  the  first  class 
of  integrals,  consider  the  integral  (Cm|u|Sn)  with  m  #  n. 
Integrating  by  parts. 


<Cm|u|Sn)  =  (u/CmSndu)  -  I  /C^du'  du 


(C.l) 


Reid  and  Harris  list 


<Am4-yn4>  'cmSndu  =  Xm4/CmSndu 


4  W  W  M  - 

“  v*n  f  ^m®n^u  f  CmSndu— /  CmSndu 


(C.  2) 


The  indefinite  integrals  on  the  right  hand  side  of  (C.2) 
are  also  evaluated  by  Reid  and  Harris.  Inserting  into 
(C.2)  and  then  (C.l),  and  remembering  that 


Cm(±|)  =  Sm(±|)  =  Cm(±^)  =  Sjfi(±|)  =  0  (C.  3) 
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we  have  the  desired  result 


««MSn>  -  «*>$  -  Z4S" 


urr 


.  4<Xm4  +  wn4)  Cm 
+  - <\m*  '-'yn4^ 


(C.  4) 


where  Cm  and  Sn  are  evaluated  at  u  «  j.  Completing  the 
other  integrals  from  the  first  class  in  a  similar  manner 
gives  (for  m  ^  n) 


•  _  *> 


2  Am4 
2 

*m4~*n ' 
Sm  Sm  , 


!i*m 


PnT-Pn- 


,c»hs»i  -  uffliv 


+  Xm2  tanh2 

(C.  5) 

(cm  Cn  - 

c;  Cn) 

(C.  6) 

Pm2  coth2 

(C.7) 

.  H  Hi 

(sm  sn  ~ 

Hi  a 

sm  Sn) 

(C.  8) 

f  m 

(C.  9) 

a  at 

5  Cm  Cm 

5  +.nnh2  *m 

(C.10) 

4Xma 

2X^7  tanh  2 

w  /// 

8Cm  Cn 


-  (Cm|u2|Cn)  =  ■{-Xin4-Xn'4)  '£ 

+  (Xm4-Xn4)3  MlOXm4+6Xn4)Cm  Cn  +  (6Xm4+10Xn4) 

(C.ll) 


♦/  ♦ 
Cm  Cn} 


(SmluJISp)  -  ij  *  ^M-fm  .  coth2-mz 

»!  Hi 

(Smlu2lSn)  -  (^u%2 


(C. 12) 
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{-(10um4+6Mn4)s"  Sn  +  (6lJm4+10un4)Sm  S*} 


(C. 13) 


(Cm| u2  |  Cjf)  =  \  tanh2 


//  it!  ***  It 

I  9i„.,  cm  Cn  _  Cm  Cn 
(Cm  |  u2  |  Cn)  =  "2-x"3-Tn^r 


(C .  1 4 ) 


4(Xm4+Xn4)Cm  Cn 
( *m4"*n4) 2 


(C . 15 ) 


//  ry  O 

<S™|u2|Sri)  -  §  -  coth2  +  ftt-j  (C.16) 


with  the  C  and  S  functions  again  evaluated  at  u  =  ?• 

The  second  class  of  exponential  integrals  are  evaluated 
using  the  explicit  forms  of  C  and  S.  The  results  are  (with 
m  now  allowed  to  equal  n) 


(Cm  I e“au | Cn)  = 


_ 1 _  sinh  2  (*m+*n-a) 

2  cosh  cosh  -^5-  *  *m  +  *n  -  a 


+  sinh  ?(Xm~Xn  a)  +  sinh  1  (An-Am'a)  +  sinh  |(Xm+Xn+a) 


” *n 


xn  -Am  a 


+  *n  +a 


_ 1 _  2  ( Am_a)  cos— jsinh  +  2  Ansin-^coshXm^-a 

2  cosh  Xjfl  cos  2lu  1 - :: - 

22  *n  +  ( *m“a) 2 

2  ( Am+a) cos-^sinh-^mga  +  2\n  sin  cosh 

^  +  ( *m+a ) 2 

^  (*m_*n)sin  Xlfl2  cosh  j  +  cos  sinh  y 


a2+  (Am+\n)^ 
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(Sml e_au | Sn)  = 


2  sinh  sinh  ™  sinh  j  (Pm+Pn-a) 


Pm+Pn~a 


_  sinh  j  (pm-pn-a)  _  sinh  j  (Pn-Pm~a)  +  sinh  j  (Pm+Pn+a) 


Pm-Pn-a 


Pn-Pm-a 


Pm+Pn+a 


2  sinh  sin 


Pn  {' 


2(Pm~a)sin  cosh  -  2pncos^sinh-^:!Jfl^- 


Pn2  +  (Pm"3)' 


+  2(ym+a)sin  cosh  -Ja-^a  -  2yn  cos  sinh 

Un2  +  ( Pm+a ) ^  ^ 

l  (Pm"Pn)  sin  cosh  ^  +  a  cos  -Vffl+Vn  sinh  ^ 


sin-^sin-^j 


a2  +  (Pm“Pn)2 


bm+Mn)  sin  HmtHn  cosh  j  +  a  cos  sinh  j 

a^  +  Twin+ynT^  ^ 

^  2(yn-a)sin  ^  cosh  -  2ym  cos-^sinh^11^ 


2  sinh  sin 


Pm 


^  +  (pn~a) 2 


2(pn+a)  sin-^S  cosh  -  nta  -  2pm  cos  sinh  |JaXa 

- f - - — - — j — — - ± - — }  (C.  18) 

Pm^  +  ( Pn+a) ^ 

!  sinh  ^  (xm  +  un-a) 


(Cm|e  |sn)  2  s.nh  ^  cosh  Xj  '  *m+Pn-a 

sinh  ^  (pn-Anra)  _  sinh  j  (xm-pn-a)  _  sinh  jUni+Pn+a) 


Pn-^m~a 


*m-Pn“a 

X 


^m+Pn+a 


2(yn-a)cos  sinh—^  +  2Xmsin  cosh 


2  cos  — ^  sinh 


*m2  +  (vn-a)2 
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_  2(un+a)  cos  sinh  An+l  +  2Xm  sin  cosh  -Hn+a 


Am2  +  (pn+a)2 


2  ( Xm-a)  sin-^&cosh-^3  -  2pnCos^4^sinhAilX 


Xm_a 


2  cosh  — sin 


^n2  +  <xm-a)2 


_  2(Xm+a)  sin  ^  cosh  -  2gra  cos  ^  sinh  — ?  } 


Xjn+a 


Pn2  "*■  (  Xm+a )  2 


cos  sin  An. 


1  ~  a  sin  Mn^mcosh  j  +  ( pn-Xm)  cos-^y^sinh 


2  +  ^  Pn'^ra) 2 


+  ~  a  sin  im±im  cosh  |  +  (pn+Xm)  cos  An+Am  sinh  J 


(C . 19 ) 


a2  +  (pn+*m)2 


(Cm|e"au|Cn)  = 


Amxn 


-t —  sinh  ■=■  (Xm+Xn  —  a) 

A  ii  i  ^ 


2  cosh  -w  cosh  — 4  {-  , 

a  2  *m  +  *n  _a 


_  sinh  ^(Xm-Xn-a)  sinh  |-(Xn-Xm-a)  sinh  ^  (*m+*n+a) 

xm-*n-a  xn_xm~a  +  Xm+x„+a  ^ 


'nr  *n 


- - -  2  ( Xm-a)  sin  ^cosh  2Xncos-^  sinh  AnCS. 

2cosh-^»cos^  { - ± - ± _ ” _ _ 2__ 

*n2  +  ( Xm-a) 2 


+  2  (Xm+a)  sin  An  coah  Amjl  -  2xn  cos  sinh  AmJ 


Xm+a 


Xn2  +  (Xm+a)2 


1 


-Ami o_ 

.v  An- 


2(Xn-a)sin  cosh  -■nT.ff  -  2 Xmcos  sinh 


(8kN 


+  2  ( An+a)  sin  cosh  ~  2xm  cos  ^  sinh  ^±1 

*m2  +  ( ^n+a) 2 


+  - -x-£~ - r—  (xm-^n)sin±m/'A~n  cosh  f  + 

cos  cos  { - - - 1 _ - _ 


•  —  ^  m  “  ^  ? 


^  +  a  cos  --ULliLll  sinh 


a2  +  ( 2 
_  s^n  cosh  §■  +  a  cos  *n  sinh  ^ 


(C.20) 


a2  +  Um+*n)2 


(Sm|e-au|Sn)  = 


MmWn 


2  sinh  sinh 


sinh  =  (Pm+vn~a) 

{ - - — — - 

Mm  +  t'n-9 


+  sinh  2  ^m_lJn_a^  +  sinh  ^  (Mn-Mm-a)  +  sinh  ^  (un  +Pm+a  ) 


Mm~Mn“a 


Mn-Mm-a 


Mm+Mn+a 


^—-0 — - -  2(um-a)  cos  ^  sinh  +  2pnsin-^cosh-^ffl- 


-a 


2  sinh  ^  sin  ^  { 


Mn2  +  (um-a)2 


+  2(un+a)  cos  sinh  ~n^a+  2pm  sin  cosh  pn+a 

wm2  +  (un+a>2 


- —  P-[? -  2  ( un-a )  cos-^sinh^"8  +  2  wmsin-^^cosh-^C^ 

2  sinh  ^  sin  ^  2  m - —L 

2  2  unr  +  (yn~a>2 


+  2(pn+a)  cos  sinh  Pnta+  2pm  sin  cosh  -J1y^ 


Mm 


2  +  (ijn+a)2 


Wn 


“nTV 


4  U 

sin  -4j  sin  -jy 


(4m~4n)sin  “‘g  “cosh  |  +  a  cos-2^— ^sinh  | 


2  A  ,  x 2 

a  +  (lx  -a  ) 
v^m  Kn' 


66 


f 


li  +  (1  U  +jl 

+  (^m+Un)  sin  — ^  cosh  |  +  a  cos  mg  ~  sinh  | 


2  ,  ,2 
a  +  (|i  +(i  ) 
m  n 


and 


tc ' !  p*au  I  o  -  \  -  _ kfnxn _  sinh  (wn+um-a) 

'cml e  I  bn'  2  cosh  M  sinh  Ah  { - - - 

2  2  ^n+^m-3 


+  sinh  ^(Am-nn-a)  _  sinh  -j  (Mn“Am-a)  _  sinh  ^(un+Am+a) 


Am-un-a 


Mn- Am-a 


Wn+ Am+a 


Am^n _ 

2sin-u^cosh-^ 


2(Am-a)cos  •^&sinh-^m^  +  2pnsin-^?cosh 


Jn2  +  ( Am-a) 2 


2(Am+a)  cos  ^  sinh  Am  ■--+  2yn  sin  cosh 


_  j  _  U n _ %_  Am+& 

-2-) 


v*n 


2  +  ( Am+a) 2 


2  ( pn~a)  sin-^cosh-^11^  -  2  Amcos-^Ssinh 


_ AmHn _  , 

2cos^sinhA^D.  1 


-a 


Am2  +  (un*a)2 


2(yn+a)  sin  ^  cosh  Pm^—  -  2Am  cos  sinh  Ufl^a 


Am2  +  ( Pn+a) 2 


1 


— -  _  a  sin  ■  mI-n  cosh  §^+(Am-yn)  cos  All^n  sinh 

osAHlsiniiS  (  - 

2  ^  a2  +  ( Am_un) 2 


f  ^m-rWn; 


+  -  a  sin  Xm^-n  cosh  j 

a2  +  (Am+wn)2 


All  integrals  listed  in  this  appendix  have  been 
verified  using  a  Simpson  approximation  method  for  (m,n) 
equal  to  (1,1),  (2,2),  and  (1,2). 
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